
plot.sim = function(sim, y.range.spud=NULL, cex.leg=0.85) {
  par(mfrow=c(2,2))
  
  plot(wti ~ date, data=sim$PricePaths[[1]], lty=1, col='black', lwd=2, type='l',
       xlab='Date',ylab='$/barrel of oil equivalent', main='Oil and Gas Prices', ylim=c(0,80))
  lines(I(hh*6) ~ date, data=sim$PricePaths[[1]], lty=1, col='blue', lwd=2)
  legend('top', legend=c('WTI','Henry Hub'),
         lty=1, lwd=2, col=c('black','blue'), cex=cex.leg, bty='n',
         y.intersp=0.5)
  grid(nx=NA, ny=NULL)
  
  plot(Oil_Nonfederal ~ date, data=sim$spuds, lty=1, col=my.cols[1], lwd=2, type='l', 
       ylim=y.range.spud, xlab='Date', ylab='Wells drilled per month', main='Wells Drilled')
  grid(nx=NA, ny=NULL)
  lines(Oil_Federal ~ date, data=sim$spuds, lty=1, col=my.cols[2], lwd=2)
  lines(Gas_Nonfederal ~ date, data=sim$spuds, lty=1, col=my.cols[3], lwd=2)
  lines(Gas_Federal ~ date, data=sim$spuds, lty=1, col=my.cols[4], lwd=2)
  legend('top',legend=c('Oil Wells, Nonfederal', 'Oil Wells, Federal',
                        'Gas Wells, Nonfederal', 'Gas Wells, Federal'),
         ncol=2, cex=cex.leg, bty='n', x.intersp=1, y.intersp=0.5, text.width=2100,
         col=my.cols, lwd=3, lty=1)
  
  plot(liq_Total ~ date, data=sim$OilProduction, type='l', lwd=2, main='Oil Production', 
       xlab='Date', ylab='mb/d')
  grid(nx=NA, ny=NULL)
  legend('right', legend='Total', col='black', lwd=2, lty=1, bty='n')
  lines(Oil_Nonfederal ~ date, data=sim$OilProduction, lty=1, col=my.cols[1], lwd=2)
  lines(Oil_Federal ~ date, data=sim$OilProduction, lty=1, col=my.cols[2], lwd=2)
  lines(Gas_Nonfederal ~ date, data=sim$OilProduction, lty=1, col=my.cols[3], lwd=2)
  lines(Gas_Federal ~ date, data=sim$OilProduction, lty=1, col=my.cols[4], lwd=2)
  
  plot(gas_Total ~ date, data=sim$GasProduction, type='l', lwd=2, main='Gas Production',
       xlab='Date', ylab='bcf/d')
  grid(nx=NA, ny=NULL)
  lines(Oil_Nonfederal ~ date, data=sim$GasProduction, lty=1, col=my.cols[1], lwd=2)
  lines(Oil_Federal ~ date, data=sim$GasProduction, lty=1, col=my.cols[2], lwd=2)
  lines(Gas_Nonfederal ~ date, data=sim$GasProduction, lty=1, col=my.cols[3], lwd=2)
  lines(Gas_Federal ~ date, data=sim$GasProduction, lty=1, col=my.cols[4], lwd=2)
}

plot.us.sim = function(sim, price.range=c(0,120), y.range.spud=c(0,2400), date.range=NULL,
                       y.range.oil=c(0,20), y.range.gas=c(0,150), cex.leg=0.95) {
  if (is.null(date.range)) date.range = sim$spuds[,date]
  par(mfrow=c(3,4))
  
  plot(wti ~ date, data=sim$PricePaths[[1]][date %in% date.range], lty=1, col='black', lwd=2, type='l',
       xlab='',ylab='$/barrel (left) or $/mmbtu (right)', main='Oil and Gas Prices', ylim=price.range)
  grid(nx=NA, ny=NULL)
  abline(v=seq.Date(as.Date('2020-01-01'),as.Date('2050-01-01'), by='5 years'), col='lightgray', lty=3)
  abline(v=sim.start.date, lty=1, col='lightgray')
  lines(wti ~ date, data=sim$PricePaths[[1]][date %in% date.range], lty=1, col='black', lwd=2)
  par(new=TRUE)
  plot(hh ~ date, data=sim$PricePaths[[1]][date %in% date.range], lty=1, col='blue', 
       lwd=2, type='l', yaxt='n', xaxt='n', ylim=price.range/20,
       ylab='', xlab='')
  axis(4)
  legend('top', legend=c('WTI (left)','Henry Hub (right)'),
         lty=1, lwd=2, col=c('black','blue'), cex=cex.leg, bty='n',
         y.intersp=0.75)
  
  # Wells Drilled, Onshore
  plot(Oil_Nonfederal_Onshore ~ date, data=sim$spuds[date %in% date.range], lty=1, col=my.cols[1], lwd=2, type='l', 
       ylim=y.range.spud, xlab='', ylab='Wells drilled per month', main='Onshore Wells Drilled')
  abline(v=seq.Date(as.Date('2010-01-01'),as.Date('2050-01-01'), by='5 years'), col='lightgray', lty=3)
  grid(nx=NA, ny=NULL)
  lines(Oil_Nonfederal_Onshore ~ date, data=sim$spuds[date %in% date.range], lty=1, col=my.cols[1], lwd=2)
  lines(Oil_Federal_Onshore ~ date, data=sim$spuds[date %in% date.range], lty=1, col=my.cols[2], lwd=2)
  lines(Gas_Nonfederal_Onshore ~ date, data=sim$spuds[date %in% date.range], lty=1, col=my.cols[3], lwd=2)
  lines(Gas_Federal_Onshore ~ date, data=sim$spuds[date %in% date.range], lty=1, col=my.cols[4], lwd=2)
  abline(v=sim.start.date, lty=1, col='lightgray')
  
  # Wells Drilled, Offshore
  plot(Oil_Federal_Offshore ~ date, data=sim$spuds[date %in% date.range], lty=1, col=my.cols[2], lwd=2, main='Offshore Wells Drilled', ylab='Wells Drilled per month', type='l')
  abline(v=seq.Date(as.Date('2020-01-01'),as.Date('2050-01-01'), by='5 years'), col='lightgray', lty=3)
  grid(nx=NA, ny=NULL)
  lines(Oil_Nonfederal_Offshore ~ date, data=sim$spuds[date %in% date.range], lty=1, col=my.cols[1], lwd=2)
  lines(Oil_Federal_Offshore ~ date, data=sim$spuds[date %in% date.range], lty=1, col=my.cols[2], lwd=2)
  lines(Gas_Nonfederal_Offshore ~ date, data=sim$spuds[date %in% date.range], lty=1, col=my.cols[3], lwd=2)
  lines(Gas_Federal_Offshore ~ date, data=sim$spuds[date %in% date.range], lty=1, col=my.cols[4], lwd=2)
  abline(v=sim.start.date, lty=1, col='lightgray')
  
  # Blank panel
  plot.new()
  legend('left',legend=c('Oil Wells, Nonfederal', 'Oil Wells, Federal',
                         'Gas Wells, Nonfederal', 'Gas Wells, Federal'),
         ncol=1, cex=cex.leg*1.5, bty='n', x.intersp=1, y.intersp=1, text.width=0.5*strwidth('Oil Wells, Nonfederal'),
         col=my.cols, lwd=3, lty=1)
  
  # Oil Production, New Wells
  plot(liq_Total ~ date, data=sim$OilProduction_NewWells[date %in% date.range], type='l', lwd=2, main='Oil Production\nNew Wells', 
       xlab='', ylab='mb/d', ylim=y.range.oil)
  abline(v=seq.Date(as.Date('2010-01-01'),as.Date('2050-01-01'), by='5 years'), col='lightgray', lty=3)
  grid(nx=NA, ny=NULL)
  legend('top', legend='Total', col='black', lwd=2, lty=1, bty='n')
  lines(liq_Total ~ date, data=sim$OilProduction_NewWells[date %in% date.range], lty=1, col='black', lwd=2)
  lines(Oil_Nonfederal_Onshore ~ date, data=sim$OilProduction_NewWells[date %in% date.range], lty=1, col=my.cols[1], lwd=2)
  lines(Oil_Federal_Onshore ~ date, data=sim$OilProduction_NewWells[date %in% date.range], lty=1, col=my.cols[2], lwd=2)
  lines(Gas_Nonfederal_Onshore ~ date, data=sim$OilProduction_NewWells[date %in% date.range], lty=1, col=my.cols[3], lwd=2)
  lines(Gas_Federal_Onshore ~ date, data=sim$OilProduction_NewWells[date %in% date.range], lty=1, col=my.cols[4], lwd=2)
  
  lines(Oil_Nonfederal_Offshore ~ date, data=sim$OilProduction_NewWells[date %in% date.range], lty=2, col=my.cols[1], lwd=2)
  lines(Oil_Federal_Offshore ~ date, data=sim$OilProduction_NewWells[date %in% date.range], lty=2, col=my.cols[2], lwd=2)
  lines(Gas_Nonfederal_Offshore ~ date, data=sim$OilProduction_NewWells[date %in% date.range], lty=2, col=my.cols[3], lwd=2)
  lines(Gas_Federal_Offshore ~ date, data=sim$OilProduction_NewWells[date %in% date.range], lty=2, col=my.cols[4], lwd=2)
  abline(v=sim.start.date, lty=1, col='lightgray')
  
  # Oil Production, Existing Wells
  plot(liq_Total ~ date, data=sim$OilProduction_ExWells[date %in% date.range], type='l', lwd=2, main='Oil Production\nExisting Wells', 
       xlab='', ylab='mb/d', ylim=y.range.oil)
  abline(v=seq.Date(as.Date('2010-01-01'),as.Date('2050-01-01'), by='5 years'), col='lightgray', lty=3)
  grid(nx=NA, ny=NULL)
  lines(liq_Total ~ date, data=sim$OilProduction_ExWells[date %in% date.range], lty=1, col='black', lwd=2)
  lines(Oil_Nonfederal_Onshore ~ date, data=sim$OilProduction_ExWells[date %in% date.range], lty=1, col=my.cols[1], lwd=2)
  lines(Oil_Federal_Onshore ~ date, data=sim$OilProduction_ExWells[date %in% date.range], lty=1, col=my.cols[2], lwd=2)
  lines(Gas_Nonfederal_Onshore ~ date, data=sim$OilProduction_ExWells[date %in% date.range], lty=1, col=my.cols[3], lwd=2)
  lines(Gas_Federal_Onshore ~ date, data=sim$OilProduction_ExWells[date %in% date.range], lty=1, col=my.cols[4], lwd=2)
  
  lines(Oil_Nonfederal_Offshore ~ date, data=sim$OilProduction_ExWells[date %in% date.range], lty=2, col=my.cols[1], lwd=2)
  lines(Oil_Federal_Offshore ~ date, data=sim$OilProduction_ExWells[date %in% date.range], lty=2, col=my.cols[2], lwd=2)
  lines(Gas_Nonfederal_Offshore ~ date, data=sim$OilProduction_ExWells[date %in% date.range], lty=2, col=my.cols[3], lwd=2)
  lines(Gas_Federal_Offshore ~ date, data=sim$OilProduction_ExWells[date %in% date.range], lty=2, col=my.cols[4], lwd=2)
  abline(v=sim.start.date, lty=1, col='lightgray')
  
  # Total liq, new+existing, by type
  plot(liq_Total ~ date, data=sim$OilProduction_AllWells[date %in% date.range], type='l', lwd=2, main='Oil Production\nAll Wells', 
       xlab='', ylab='mb/d', ylim=y.range.oil)
  abline(v=seq.Date(as.Date('2010-01-01'),as.Date('2050-01-01'), by='5 years'), col='lightgray', lty=3)
  grid(nx=NA, ny=NULL)
  lines(liq_Total ~ date, data=sim$OilProduction_AllWells[date %in% date.range], lty=1, col='black', lwd=2)
  lines(Oil_Nonfederal_Onshore ~ date, data=sim$OilProduction_AllWells[date %in% date.range], lty=1, col=my.cols[1], lwd=2)
  lines(Oil_Federal_Onshore ~ date, data=sim$OilProduction_AllWells[date %in% date.range], lty=1, col=my.cols[2], lwd=2)
  lines(Gas_Nonfederal_Onshore ~ date, data=sim$OilProduction_AllWells[date %in% date.range], lty=1, col=my.cols[3], lwd=2)
  lines(Gas_Federal_Onshore ~ date, data=sim$OilProduction_AllWells[date %in% date.range], lty=1, col=my.cols[4], lwd=2)
  
  lines(Oil_Nonfederal_Offshore ~ date, data=sim$OilProduction_AllWells[date %in% date.range], lty=2, col=my.cols[1], lwd=2)
  lines(Oil_Federal_Offshore ~ date, data=sim$OilProduction_AllWells[date %in% date.range], lty=2, col=my.cols[2], lwd=2)
  lines(Gas_Nonfederal_Offshore ~ date, data=sim$OilProduction_AllWells[date %in% date.range], lty=2, col=my.cols[3], lwd=2)
  lines(Gas_Federal_Offshore ~ date, data=sim$OilProduction_AllWells[date %in% date.range], lty=2, col=my.cols[4], lwd=2)
  abline(v=sim.start.date, lty=1, col='lightgray')
  
  # Blank panel
  plot.new()
  legend('left',legend=c('Onshore','Offshore'),
         ncol=1, cex=cex.leg*1.75, bty='n', x.intersp=1, y.intersp=1,
         text.width=0.75*strwidth('Onshore'),
         col=c('black'), lwd=2, lty=1:2)
  
  # Gas Production, New Wells
  plot(gas_Total ~ date, data=sim$GasProduction_NewWells[date %in% date.range], type='l', lwd=2, main='Gas Production\nNew Wells',
       xlab='', ylab='bcf/d', ylim=y.range.gas)
  abline(v=seq.Date(as.Date('2010-01-01'),as.Date('2050-01-01'), by='5 years'), col='lightgray', lty=3)
  grid(nx=NA, ny=NULL)
  legend('top', legend='Total', col='black', lwd=2, lty=1, bty='n')
  lines(gas_Total ~ date, data=sim$GasProduction_NewWells[date %in% date.range], lty=1, col='black', lwd=2)
  lines(Oil_Nonfederal_Onshore ~ date, data=sim$GasProduction_NewWells[date %in% date.range], lty=1, col=my.cols[1], lwd=2)
  lines(Oil_Federal_Onshore ~ date, data=sim$GasProduction_NewWells[date %in% date.range], lty=1, col=my.cols[2], lwd=2)
  lines(Gas_Nonfederal_Onshore ~ date, data=sim$GasProduction_NewWells[date %in% date.range], lty=1, col=my.cols[3], lwd=2)
  lines(Gas_Federal_Onshore ~ date, data=sim$GasProduction_NewWells[date %in% date.range], lty=1, col=my.cols[4], lwd=2)
  
  lines(Oil_Nonfederal_Offshore ~ date, data=sim$GasProduction_NewWells[date %in% date.range], lty=2, col=my.cols[1], lwd=2)
  lines(Oil_Federal_Offshore ~ date, data=sim$GasProduction_NewWells[date %in% date.range], lty=2, col=my.cols[2], lwd=2)
  lines(Gas_Nonfederal_Offshore ~ date, data=sim$GasProduction_NewWells[date %in% date.range], lty=2, col=my.cols[3], lwd=2)
  lines(Gas_Federal_Offshore ~ date, data=sim$GasProduction_NewWells[date %in% date.range], lty=2, col=my.cols[4], lwd=2)
  abline(v=sim.start.date, lty=1, col='lightgray')
  
  # Gas Production, Existing Wells
  plot(gas_Total ~ date, data=sim$GasProduction_ExWells[date %in% date.range], type='l', lwd=2, main='Gas Production\nExisting Wells',
       xlab='', ylab='bcf/d', ylim=y.range.gas)
  abline(v=seq.Date(as.Date('2010-01-01'),as.Date('2050-01-01'), by='5 years'), col='lightgray', lty=3)
  grid(nx=NA, ny=NULL)
  lines(gas_Total ~ date, data=sim$GasProduction_ExWells[date %in% date.range], lty=1, col='black', lwd=2)
  lines(Oil_Nonfederal_Onshore ~ date, data=sim$GasProduction_ExWells[date %in% date.range], lty=1, col=my.cols[1], lwd=2)
  lines(Oil_Federal_Onshore ~ date, data=sim$GasProduction_ExWells[date %in% date.range], lty=1, col=my.cols[2], lwd=2)
  lines(Gas_Nonfederal_Onshore ~ date, data=sim$GasProduction_ExWells[date %in% date.range], lty=1, col=my.cols[3], lwd=2)
  lines(Gas_Federal_Onshore ~ date, data=sim$GasProduction_ExWells[date %in% date.range], lty=1, col=my.cols[4], lwd=2)
  
  lines(Oil_Nonfederal_Offshore ~ date, data=sim$GasProduction_ExWells[date %in% date.range], lty=2, col=my.cols[1], lwd=2)
  lines(Oil_Federal_Offshore ~ date, data=sim$GasProduction_ExWells[date %in% date.range], lty=2, col=my.cols[2], lwd=2)
  lines(Gas_Nonfederal_Offshore ~ date, data=sim$GasProduction_ExWells[date %in% date.range], lty=2, col=my.cols[3], lwd=2)
  lines(Gas_Federal_Offshore ~ date, data=sim$GasProduction_ExWells[date %in% date.range], lty=2, col=my.cols[4], lwd=2)
  abline(v=sim.start.date, lty=1, col='lightgray')
  
  # Total gas, new+existing, by type
  plot(gas_Total ~ date, data=sim$GasProduction_AllWells[date %in% date.range], type='l', lwd=2, main='Gas Production\nAll Wells', 
       xlab='', ylab='bcf/d', ylim=y.range.gas)
  abline(v=seq.Date(as.Date('2010-01-01'),as.Date('2050-01-01'), by='5 years'), col='lightgray', lty=3)
  grid(nx=NA, ny=NULL)
  lines(gas_Total ~ date, data=sim$GasProduction_AllWells[date %in% date.range], lty=1, col='black', lwd=2)
  lines(Oil_Nonfederal_Onshore ~ date, data=sim$GasProduction_AllWells[date %in% date.range], lty=1, col=my.cols[1], lwd=2)
  lines(Oil_Federal_Onshore ~ date, data=sim$GasProduction_AllWells[date %in% date.range], lty=1, col=my.cols[2], lwd=2)
  lines(Gas_Nonfederal_Onshore ~ date, data=sim$GasProduction_AllWells[date %in% date.range], lty=1, col=my.cols[3], lwd=2)
  lines(Gas_Federal_Onshore ~ date, data=sim$GasProduction_AllWells[date %in% date.range], lty=1, col=my.cols[4], lwd=2)
  
  lines(Oil_Nonfederal_Offshore ~ date, data=sim$GasProduction_AllWells[date %in% date.range], lty=2, col=my.cols[1], lwd=2)
  lines(Oil_Federal_Offshore ~ date, data=sim$GasProduction_AllWells[date %in% date.range], lty=2, col=my.cols[2], lwd=2)
  lines(Gas_Nonfederal_Offshore ~ date, data=sim$GasProduction_AllWells[date %in% date.range], lty=2, col=my.cols[3], lwd=2)
  lines(Gas_Federal_Offshore ~ date, data=sim$GasProduction_AllWells[date %in% date.range], lty=2, col=my.cols[4], lwd=2)
  abline(v=sim.start.date, lty=1, col='lightgray')
  
  # Focus on offshore, all wells.
  plot(I(Oil_Federal_Offshore) ~ date, data=sim$OilProduction_AllWells[date %in% date.range], type='l',
       lty=1, col=my.cols[2], lwd=2, main='Offshore Oil & Gas  Production\nAll Wells', ylab='mb/d (left) or bcf/d (right)', ylim=c(0,3))
  abline(v=seq.Date(as.Date('2010-01-01'),as.Date('2050-01-01'), by='5 years'), col='lightgray', lty=3)
  grid(nx=NA, ny=NULL)
  lines(I(Oil_Nonfederal_Offshore) ~ date, data=sim$OilProduction_AllWells[date %in% date.range], lty=1, col=my.cols[1], lwd=2)
  lines(I(Oil_Federal_Offshore) ~ date, data=sim$OilProduction_AllWells[date %in% date.range], lty=1, col=my.cols[2], lwd=2)
  lines(I(Gas_Nonfederal_Offshore) ~ date, data=sim$OilProduction_AllWells[date %in% date.range], lty=1, col=my.cols[3], lwd=2)
  lines(I(Gas_Federal_Offshore) ~ date, data=sim$OilProduction_AllWells[date %in% date.range], lty=1, col=my.cols[4], lwd=2)
  lines(I(Oil_Federal_Offshore+Gas_Federal_Offshore) ~ date, data=sim$OilProduction_AllWells[date %in% date.range], lty=1, col='black', lwd=2)
  par(new=TRUE)
  plot(I(Oil_Nonfederal_Offshore) ~ date, data=sim$GasProduction_AllWells[date %in% date.range], type='l', lty=3, col=my.cols[1], lwd=2,
       xaxt='n',yaxt='n',ylab='',xlab='', ylim=c(0,8))
  axis(4)
  lines(I(Oil_Federal_Offshore) ~ date, data=sim$GasProduction_AllWells[date %in% date.range], lty=2, col=my.cols[2], lwd=2)
  lines(I(Gas_Nonfederal_Offshore) ~ date, data=sim$GasProduction_AllWells[date %in% date.range], lty=2, col=my.cols[3], lwd=2)
  lines(I(Gas_Federal_Offshore) ~ date, data=sim$GasProduction_AllWells[date %in% date.range], lty=2, col=my.cols[4], lwd=2)
  lines(I(Oil_Federal_Offshore+Gas_Federal_Offshore) ~ date, data=sim$GasProduction_AllWells[date %in% date.range], lty=2, col='black', lwd=2)
  legend('top', legend=c('Oil Production (left)','Gas Production (right)'), lty=1:2, lwd=2, col='black',
         y.intersp = 1, bty='n', horiz=FALSE, text.width = 0.95*strwidth('Oil Production (left)'),
         cex=cex.leg)
  abline(v=sim.start.date, lty=1, col='lightgray')
}



plot.global.sim = function(sim, price.range=NULL, date.range=NULL, cex.leg=1, plot.inven=FALSE) {
  if (is.null(date.range)) date.range = sim$global.sim[, unique(date)]
  if (is.null(price.range)) price.range = sim$global.sim[, c(0, 1.1*max(wti))]
  setkey(sim$global.sim, date)
  
  date.grid = seq.Date(as.Date('2020-01-01'), as.Date('2050-01-01'), by='5 years')
  
  # Plot prices
  par(mfrow=ifelse(rep(plot.inven, 2), c(1,4), c(1,3)))
  plot(brent ~ date, data=sim$global.sim[J(date.range)], lty=1, col='black', lwd=2, type='l',
       xlab='',ylab='$/barrel (left) or $/mmbtu (right)', main='Oil and Gas Prices', cex.main=1.5, ylim=price.range,
       cex.axis=1.5, cex.lab=1.5)
  lines(brent.base ~ date, data=sim$global.sim[J(date.range)], col='black', lwd=2, lty=2)
  grid(nx=NA, ny=NULL); abline(v=date.grid, lty=3, col='lightgray')
  abline(h=0); abline(v=sim.start.date, col='lightgray')
  
  par(new=TRUE)
  plot(hh ~ date, data=sim$global.sim[J(date.range)], lty=1, col='blue', 
       lwd=2, type='l', yaxt='n', xaxt='n', ylim=price.range/20,
       ylab='', xlab='')
  
  lines(hh.base ~ date , data=sim$global.sim[J(date.range)], col='blue', lwd=2, lty=2)
  axis(4, cex.axis=1.5)
  legend('top', legend=c('Brent, $/bbl (left)','Henry Hub, $/mmbtu (right)'),
         lty=1, lwd=2, col=c('black','blue'), cex=cex.leg, bty='n',
         horiz=F)
  
  
  # Plot US/ROW/Total Oil production
  plot(Tot.Oil.Supply ~ date, data=sim$global.sim[J(date.range)], type='l', ylim=c(-10,150), 
       lwd=2, col='black', lty=1, main='Crude Oil Supply & Demand', cex.main=1.5,
       ylab='mb/d', xlab='',
       cex.axis=1.5, cex.lab=1.5)
  
  lines(Tot.Oil.Demand ~ date, data=sim$global.sim[J(date.range)], lty=2, lwd=2, col='green')
  lines(US.oil.supply ~ date, data=sim$global.sim[J(date.range)], lty=1, lwd=2, col='blue')
  lines(row.oil.supply ~ date, data=sim$global.sim[J(date.range)], lty=1, lwd=2, col='red')
  lines(Excess.Oil.Demand ~ date, data=sim$global.sim[J(date.range)], lty=1, lwd=2, col='purple')
  legend('top', legend=c('Total Demand','Total Supply','US Supply','ROW Supply'), horiz=F, 
         lty=c(2,1,1,1), lwd=2, cex=cex.leg, col=c('green','black','blue','red'), bty='n')
  grid(nx=NA, ny=NULL); abline(v=date.grid, lty=3, col='lightgray')
  abline(h=0); abline(v=sim.start.date, col='lightgray')
  abline(h=c(10,20), lty=3, col='lightgray'); axis(2, at=c(10,20), cex.axis=1.5)
  
  # Plot US/ROW/Total Gas production
  plot(Tot.Gas.Supply ~ date, data=sim$global.sim[J(date.range)], type='l', ylim=c(-40,600),
       lwd=2, col='black', lty=1, main='Gas Supply & Demand', cex.main=1.5,
       ylab='bcf/d', xlab='',
       cex.axis=1.5, cex.lab=1.5)
  
  lines(Tot.Gas.Demand ~ date, data=sim$global.sim[J(date.range)], lty=2, lwd=2, col='green')
  lines(US.gas.supply ~ date, data=sim$global.sim[J(date.range)], lty=1, lwd=2, col='blue')
  lines(row.gas.supply ~ date, data=sim$global.sim[J(date.range)], lty=1, lwd=2, col='red')
  lines(Excess.Gas.Demand ~ date, data=sim$global.sim[J(date.range)], lty=1, lwd=2, col='purple')
  legend('top', legend=c('Total Demand','Total Supply','US Supply','ROW Supply'), horiz=F, 
         lty=c(2,1,1,1), lwd=2, cex=cex.leg, col=c('green','black','blue','red'), bty='n')
  grid(nx=NA, ny=NULL); abline(v=date.grid, lty=3, col='lightgray')
  abline(h=0); abline(v=sim.start.date, col='lightgray')
  
  # Plot inventories
  if (plot.inven) {
    plot(Oil.Inventories ~ date, data=sim$global.sim[J(date.range)], type='l', 
         ylim=c(-5000, 10000), lwd=2, col='black', lty=2, main='Inventories', ylab='Million Barrels',
         xlab='', cex.axis=1.5, cex.lab=1.5)
    par(new=TRUE)
    plot(Gas.Inventories ~ date, data=sim$global.sim[J(date.range)],
         yaxt='n', xaxt='n', ylab='', xlab='',
         type='l', ylim=c(-5000, 10000), lwd=2, col='blue', lty=2, main='', cex.main=1.5,
         cex.axis=1.5, cex.lab=1.5)
    axis(4)
    grid(); abline(h=0); abline(v=sim.start.date, col='lightgray')
  }
}
